Center for Turbulence Research 
Proceedings of the Summer Program 1992 


/ ^ 7 L99 

N9 4- 14 757 


Rapid distortion analysis and direct 
simulation of compressible homogeneous 
turbulence at finite Mach number 


By C. Cambon, 1 G. N. Coleman 2 AND N. N. Mansour 3 


The effect of rapid mean compression on compressible turbulence at a range of 
turbulent Mach numbers is investigated. Rapid distortion theory (RDT) and direct 
numerical simulation results for the case of axial (one-dimensional) compression are 
used to illustrate the existence of two distinct rapid compression regimes. These 
regimes are set by the relationships between the timescales of the mean distortion, 
the turbulence, and the speed of sound. A general RDT formulation is developed 
and is proposed as a means of improving turbulence models for compressible flows. 

1. Introduction 

This paper focuses upon the behavior of homogeneous compressible turbulence 
under the influence of rapid axial (one-dimensional) mean compression. The mo- 
tivation for this study is a need to cast light upon the physics of compressible 
turbulent flows and to improve compressible turbulence models. Our approach is 
to use both direct numerical simulations (DNS) and rapid distortion theory (RDT). 
The RDT developed in this paper is for general (those that preserve homogeneity) 
mean deformations; the resulting insight is then used to suggest improvements to 
compressible turbulence models that are applied to rapidly compressed flows. 

Earlier RDT studies of homogeneous compressible turbulence have been limited 
to either isotropic compressions (Blaisdell 1992, private communication) or the van- 
ishing turbulent Mach number limit (Durbin & Zeman 1992, hereafter referred to 
as DZ); the present investigation, therefore, attempts a more general treatment in 
that non-isotropic compressions and finite Mach numbers are considered. Some of 
our main conclusions confirm and extend those found in the recent study of shock- 
turbulence interactions by Jacquin & Cambon (1992). 

An overview of our findings follows. The RDT analysis predicts that the crucial 
parameter for turbulence subjected to rapid compression is the ratio of the mean 
deformation rate, D, to the inverse sonic timescale L/a, where L is a turbulent 
lengthscale and a is the sound speed. This parameter, DL/a , hereinafter denoted 
as Am (after DZ), is equivalent to the product of the inverse of the turbulent 
timescale, the deformation rate, and the turbulent Mach number, M t \ it defines for 
the dilatational part of the velocity field two distinct limits: the “pseudo-acoustical” 
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(nearly solenoidal) regime given by Am < 1 (that studied by DZ) and the so-called 
“pressure- released” regime with Am > 1. The term “pressure-released” is chosen 
because when Am is large, the sonic and turbulent timescales are both much larger 
than D~ l and, therefore, correlations involving the fluctuating pressure and velocity 
fields are negligible during a rapid distortion. The behavior of the solenoidal velocity 
field, according to the RDT analysis, is unaffected by the dilatational field when the 
mean flow is irrotational, and is thus independent of Am for axial compressions. Its 
history is, therefore, identical to that predicted for compression of purely solenoidal 
turbulence. In the following, we confirm these RDT predictions by comparison with 
DNS results. 

The DNS results also show that for moderate values of Am, all the double- velocity 
correlations involving the dilatational part of the turbulent velocity field remain 
weak with respect to the pure- solenoidal correlations and are in this sense similar 
to the pure solenoidal case, even for moderate compressibility. Only the Am » 1 
case is characterized by a strong amplification of the dilatational correlations. 

The moderate Am results are at first glance in conflict with recent studies of 
axially compressed turbulence (e.g., DZ, Zeman & Coleman 1992) which find un- 
expectedly large pressure-dilatation correlations in the nearly solenoidal flow. This 
led us to investigate the behavior of the pressure field, which has two roles for a 
rapid compression. On one hand, it modifies the production term in the turbulent 
kinetic energy equation by changing the Reynolds stress anisotropy through the 
classic pressure-strain rate correlation (via IIn for an axial compression in the x \ 
direction). On the other hand, the pressure is directly involved in the kinetic energy 
equation through the pressure-dilatation term, II = II;;/2. The magnitude of II u 
is found to be larger than that of II in all cases considered in this paper for a wide 
range of Mach numbers and large (but finite) compression speeds. 

Both the pressure variance and pressure-dilatation correlation from the DNS are 
found to increase with Mach number (and, therefore, with Am at a fixed mean 
distortion- to-turbulent timescale ratio) with respect to their initial values. How- 
ever, when II is compared to the production term in the turbulent kinetic energy 
transport equation, it is much smaller and has, in fact, less relative importance 
with increasing M t . This reduced relative importance of the pressure field with 
increasing compressibility is a key result of this paper and is the basis of much of 
what follows. Between the Am — > 0 and Am — > oo extremes (where the pressure- 
dilatation correlation is identically zero), II must reach a maximum; from the DNS 
results, it appears that this maximum occurs near the Am — » 0 limit at a small but 
finite value. 

In the next section, the RDT analysis is developed for compressible homogeneous 
turbulence; in §3, the theory is applied to the case of axial compression, and sep- 
arate analytic expressions for the relevant dilatational and solenoidal correlations 
for both the Am <C 1 and Am » 1 extremes are presented and compared to DNS 
results. The findings suggest that it would be appropriate for turbulence models 
to “interpolate” between the two extremes in order to accurately capture the Mt 
dependence during a rapid axial compression. We propose two methods for doing so 
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in §4 which deals with the role of, and closures for, the pressure-strain rate correla- 
tion. Section 5 considers the implications of this study for isotropically compressed 
and sheared flows, and §6 contains a recap of the main results and our conclusions. 

2. A rapid distortion analysis for compressible homogeneous turbulence 

2.1 , General considerations 

Blaisdell et al. (1991, hereafter referred to as BMR) observed that the intrin- 
sic compressibility” (the non-zero divergence) of the turbulent field often tends to 
reduce the amplification of turbulent kinetic energy produced by a mean veloc- 
ity gradient, such as a bulk compression or mean shear, with respect to the pure 
solenoidal case. This effect depends on at least three different timescales and on 
the initial turbulent field. These are the mean distortion timescale, 

To' 1 = (UijUjrf' 2 (1) 

(where Uij is the mean velocity gradient), the “turbulent decay” or “turn-over 


(where q 2 /2 is the turbulent kinetic energy and L is a lengthscale of the energy 
containing eddies), and the timescale linked to the sonic speed, 

ra" 1 = a/L. (3) 

The compression speed, r = t,/t d , is the only relevant parameter for model- 
ing homogeneous incompressible turbulence (at least for large Reynolds number). 
However, when intrinsic compressibility is considered, the ratio of the two latter 
timescales, which amounts to a turbulent Mach number M t = T a /r t , must also be 
accounted for. The magnitude of the reduction of the kinetic energy amplification 
mentioned above is, therefore, not necessarily universal, given the multi-timescale 
and initial-value nature of the problem. In fact, RDT studies of inhomogeneous 
flows even go so far as to predict an increase with M t of the kinetic energy amplifi- 
cation for turbulence under rapid (but finite) compression; these studies by Debieve 
et al (1982, hereafter referred to a s DGG) and Jacquin &: Cambon (1992) are dis- 
cussed in a following subsection, where the general RDT equations are presented 
and the reasons for the apparent growth rate versus M, discrepancy are given. 
This analysis is based on an extended Craya-Herring decomposition (Cambon 1982, 
1990; Cambon et al. 1985), which is shown to facilitate a separate investigation of 
the solenoidal and dilatational histories and provides a useful comparison to other 

approaches (e.g., BMR and DZ). l 

Some of the earlier RDT studies have apparently over-estimated the role ot the 
pressure-dilatation term, attempting to force an increased damping due to com- 
pressibility of the kinetic energy growth rate. We hope to clarify the situation here 
by separately considering various terms in one-point closure equations and thus use 
RDT as a tool for improving a model’s representation of those terms. While the 
RDT is not a model in and of itself, by improving the accuracy of crucial terms, we 
expect that it will in turn also improve the overall accuracy of the model. 
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2.2. Definitions and background 

To investigate the influence of the mean flow upon the turbulence, it is conve- 
nient to use a coordinate system x, that deforms with the mean deformation. We 
accordingly define the Lagrangian displacement tensor Fij (Eringen 1967) via 

dx ■ dx ' 

^ x i = dXj = Uidt -f FijdXj , (4) 

where £j(X, tf) is the position at time t of a fluid particle moving with the mean flow, 
which has the position Xi at the initial time t = 0. Representing the substantial 
time derivative by a superimposed dot, one has 

* dii di{ dx\ 

QXj ~~ dxi dX ~ with <Fi>(X,£ = 0,0) = Sij (5a) 

where 

< 5t > 

is the substantial derivative; we shall also have occasion to use the symbol V( ) jVt 
to denote the substantial derivative. Unless stated otherwise, the dependent vari- 
ables are assumed to be decomposed into Reynolds averaged and fluctuating com- 
ponents, as U, + t/j, where capital letters, overbars, and angle-brackets are all used 
interchangeably to denote Reynolds- (ensemble) averaged quantities, and either low- 
ercase or primed variables are used to denote fluctuating quantities. Note that F 
is a function of the stationary coordinate X, the time t, and is parameterized by 
the time (in units of t) at which the tensor is orthonormal (hence the third argu- 
ment in (5a)). For flows under mean compression, the determinant of F has special 
significance since it is equal to the volumetric ratio J. 

When the mean velocity field is irrotational, the analyses proposed (over a hun- 
dred years ago!) by Cauchy, Weber, or Kelvin for the total (mean plus fluctuating) 
vorticity can be used to give solutions for the fluctuating vorticity (w, = e ljt u*,j) 
and velocity fields: 

wj(x,f) = jF lj (X,f,0)w i (X,0) (6) 

».•(*.<) = F jr 7 I (X,<,0)ti>(X,0) + ^ ll -. (7) 

These solutions, which ultimately derive from the linearized Euler equations, remain 
approximately valid for moderately inhomogeneous flows (recall the spatial depen- 
dence of F). Eq. (6) is the classic solution of the linearized Helmholtz equation when 
the mean vorticity-fluctuating velocity term is zero (that is, for an irrotational mean 
flow). When this term is not zero, simple solutions in physical space are not possi- 
ble. Eq. (7) (also valid only for irrotational mean flows), an expression which has 
been extensively used by Goldstein (1978), contains the scalar potential <f>, which 
is directly connected to the fluctuating pressure and can be calculated once certain 
assumptions are made (e.g., that the fluctuating velocity field is solenoidal or that 
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the dilatational field is nearly acoustic). The term <f> is not the scalar potential 
arising from the Helmholtz decomposition (which we will denote <p in the following) 
because the “Ffi 1 " term in (7) contains contributions from both the solenoidal and 

dilatational velocity field. . . , 

DGG’s RDT solution of the Lagrangian transport equation for the Reynolds stress 

tensor for the case of shock wave-turbulence interaction reads 

(uiUj)(t + ) = F~l{t + ,r)(u m u n )(t-)F-\t+ ,t~), (8) 

where t~ and t + refer to positions upstream and downstream of the shock, respec- 
tively, following a mean streamline through the shock. The shock is considered as a 
pure discontinuity of the mean streamwise velocity. In other words, it is an exter- 
nal streamwise compression of infinite rate, and the associated tensor F does not 
depend on the history of the velocity gradient, but is completely characterized by 
the mean density jump or mean volumetric ratio J = Det(F), with F, } - Joabji 
through the shock. The ratio J is linked to the upstream Mach number M 0 via 

_ 2 -I- (7 - 1 )Mq ( 9) 

(7 + 1 )Mo 2 

where 7 is the ratio of specific heats, and use has been made of the classic Rankine- 
Hugoniot relations for the mean (frozen) field. A comparison of equations (7) and 
(8) shows that this approach ignores the effect of pressure (which is mediated by <j> in 
(7)); the response of the pressure fluctuations with u finite characteristic time even 
for the so-called “rapid” term is neglected compared to an infinite compression 
rate). Another idealization in the analysis of DGG, also pointed out by Lee et 
al. (1992), is that the distortion (curvature and unsteadiness) of the shock surface 
by the impinging turbulent structure is ignored. The latter issue, initially addressed 
by Ribner(1953), is not considered by the present paper. We investigate instead 
the role of the pressure field in a simpler homogeneous framework by explicitly 
defining and formalizing the range of validity of the “pressure released” regime that 
is implicit in the Debieve analysis. This paper has much in common with the recent 
analysis of the shock wave flow performed by Jacquin & Cambon (1992), in which 
the pressure-released limit was first explicitly advocated. 

Equation (7) shows that an irrotational deformation of a purely solenoidal velocity 

field is given by 

u,(x,t) = uj(x,<) = (F J 7 1 (X,t,0) U> (X,0))’ , (10) 

where to maintain = 0 we have, 

<j> ii = -{F- i \X,t,0)u } (X,0)) d (11) 

(where the s and d superscripts (and later subscripts) are understood to respectively 
refer to the solenoidal and dilatational contributions). The latter equation is an 
integral form of the Poisson equation for the fluctuating pressure, 

V 2 ^ = -(F,7 1 u i ) .. 


( 12 ) 
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For the solenoidal case, the pressure “kills off” the dilatational contribution, re- 
sulting in the lower limit of the kinetic energy growth rate caused by the mean 
compression. Conversely, in the pressure- released regime, the uf contribution in 
Eq. (11) is no longer “removed” by the pressure, producing an extra contribution 
to the solenoidal energy, which is unaffected by the dilatation field and again grows 
in accordance with Eq. (10); in other words, the compressibility leads to an increase 
in the kinetic energy growth rate. 

From this point hence, the RDT analysis will be continued under the assumption 
of flow homogeneity and make use of a spectral formalism; the Fourier wave-space 
proves to be invaluable for obtaining tractable RDT solutions. Beginning with 
Eq. (6), we shall use the Fourier space to extract the solenoidal velocity from the 
vorticity, as was done by Batchelor & Proudman (1954). Instead of solving a Poisson 
equation in physical space, we use a simple geometric wave-space projection to 
invoke the Helmholtz decomposition. 

2.3 . The mean flow 

Before turning to the turbulent fields, however, we restrict the types of mean 
deformations that are admitted by this analysis to those that preserve the homo- 
geneity of the flow. In incompressible turbulence, the constraint of maintaining 
homogeneous statistical properties leads to two conditions: the mean velocity gra- 
dient Uij must be uniform in space, and the mean flow must be a particular solution 
of the Navier- Stokes equations. The last condition amounts to an irrotational mean 
acceleration, 

vxr = o, (13) 

or that 

Uij + U i9 iUij 

is symmetric, where 

r< - (Uij + UijU^xj = Fij(t,0)Xj. (14) 

Compressibility introduces a new condition. The linearization of the momentum 
equation displays two acceleration terms. The first one is the product of mean den- 
sity and the fluctuating acceleration and leads to the same constraint mentioned 
above. The second term is the product of density fluctuation p ' by the mean acceler- 
ation T and is typically nonhomogeneous (as can be seen by the spatial dependence 
in (14)). This term can be removed, and homogeneity preserved, by neglecting the 
density fluctuation with respect to the mean density. Such an approximation (which 
is consistent with “compressed” turbulence at low Mach number) will be not used 
in this paper. Instead, we admit only mean flows without convective acceleration. 
From eqs. (5) and (14) we see that this requires 

Fij(t,0) = 6ij + Aijt , (15a) 


or 


— AtfF^ 1 — An ( Sij + Aijt) 1 . 


(15 b) 


1 1 ip i’ 


205 


RDT & DNS of compressible turbulence 

Eq. (15) is valid for an arbitrary constant (not necessarily symmetric) matrix A for 
arbitrary times, provided that the determinant of F, J , remains positive. Special 
cases of (15) have been given previously, for example, for pure strain and shear 
(BMR). A good approximation for the mean pressure P as a function of J can be 
derived from the isentropic relations. (While the isentropic relations are not strictly 
valid when M t is nonzero, DNS of finite M t turbulence under mean compression 
have shown that the deviation from the isentropic prediction is relatively small: 
for example, in Case C1DV discussed below - a rapid axial compression for initial 
Mt = 0.3 - the mean pressure at J = p(0)/p = 1/5 is within 6 percent of the 
isentropic value.) 

2.j. The fluctuating flow 

The linearized Euler equations (with p'Ti = 0) in the deforming coordinate system 


are 


Ui + Ui,jUj 



(16) 


with p = p(t) = p(0)/J(t) (recall that the dot superscript denotes a substantial 
derivative, see Eq. (5)). The linearized equations for the fluctuating pressure p and 
entropy s read (see DZ) 



= -«»,» 


i = 0, 


(17) 


where P = pRT. An investigation of the coupling between solenoidal and di- 
latational contributions to the fluctuating velocity field is conveniently done by 
transforming the variables in (the deformed coordinate) x into (three-dimensional) 
Fourier space, which we indicate either by a caret symbol or the notation U P( ).” 
The classic Helmholtz decomposition is given first in physical and then in spectral 
space as follows: 

Vi(x, t) = eijiipij + V?,.- ( is ) 

$(M) = fa - {19) 

for any vector field v. The two terms on the right-hand sides correspond to v* 
and v d , which are defined in physical space by the vector i/>, mid the ^calar po- 
tential (p. The corresponding spectral space decomposition into v' and v is given 
by the projection operators in (19), which separate the (single-component) dilata- 
tional contribution parallel to the wavevector k from the (two-component) solenoidal 
contribution in the plane normal to k. Equations (16) and (17) are easily Fourier- 
transformed; only the advection term requires particular caution: 


iii = f ( u i>t + UjjxiUij ) , 


dui 


Ui = ««,< - UijUi - 


so 
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The first and the last terms in the right-hand side of the latter equation are collec- 
tively treated as a derivative along characteristic curves, which plays the same role 
as the mean trajectories in physical space. This derivative will, therefore, also be 
represented by a superimposed dot so that 


+ Uj^ikj = 0, with solution ki = F^^OjFj. (20) 

The analogy with physical space is complete, since 

ii - Uijxj = 0, with solution Xi = Fij(t,0)Xj. (21) 

The initial k value, K, plays the same role in wave space as the Lagrangian coordi- 
nate X does in physical space. Pure kinematic distortion by advection in physical 
and spectral space are linked by a wave conservation law 

exp(iAr J a: J ) = exp^F^Xj), 

where i 2 = — 1. Accordingly, one has 

T{u x ) = Ui - UijUi ( 22 ) 


and equation (16) becomes 


Uj UljUi + — l&i ~ . 

P 


(23) 


In the latter equation, the projection operators in (19) can be used to separate 
solenoidal and dilatational contributions. We prefer to use a slightly different 
method by specifying a special frame for the solenoidal mode, according to an 
extended Craya-Herring decomposition (Cambon 1990). An orthonormal frame of 
reference (e^\e^ 2 \e^) attached to the wavevector is used with the last vector 

f o\ 

being parallel to k (e) ' = k,/k, where k is the wavevector modulus). In this local 
frame, the Fourier transform of the velocity fluctuation reads 


«i(M) = ^(M^^k) + ^ 2) (k,<)eS 2) (k) + £< 3 >(M)e[ 3) (k). (24) 

The two first terms give exactly uj, and the latter gives uf , with a minimal number 
of components and conservation of all the tensorial properties (invariants) due to the 
orthonormal properties of the local frame. Classic descriptions in terms of vorticity 
and divergence are easily recovered as 


c i),- = ifc((^^e| 2 ^ — 


(25a) 


and 


d = G til — ifc<^ 3 \ 


(256) 
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In order to remove the uncertainty regarding the azimuthal position of the solenoidal 
coordinates with respect to e< 3 \ the (e^.e^) plane is defined by choosing a fixed 
spherical coordinate polar axis n, after Herring (1974). (Craya(1958) implicitly 
used n, = S i3 and addressed only covariance matrices of the velocity field and thus 
limited the generality of his approach.) We set 

e „, = kxn md .(,) _ e .»> x ( 26 ) 

|k x n| 

Striking simplifications can be made by choosing the polar axis according to the 
symmetries of the mean flow (if any) or the statistical properties of the fluctuating 
field, retaining the full generality of the method. The equations in the local frame 
can be made nearly independent of the choice of n by using the “helical modes” 
( e (2) _ j c 0) jC ( 2 > _j_ ie| 2 ^) (which are also eigenmodes of the plane rotation matrix 
around k and of the “curl” operator) as the basis set (see Greenspan 1968, Cambon 
& Jacquin 1989, Waleffe 1992). Substituting (24) into (23) leads to the linear system 
of equations for the three components of Ui in the local Craya-Herring frame, with. 

tf a) - UirfP* + + m a3 ^ = 0 (27) 

£ (3) _ VlJ &- 3) + m 3 3 £ (3) + m 3o ^“> + i = 0. (28) 

Greek indices (indicating solenoidal space) take only the value 1 or 2, whereas the 
Latin indices range from 1 to 3 (and as in physical space, the Einstein summation 
convention is assumed). Calculation of the matrix ?n,j is straightforward, remem- 
bering to account for the rotation due to the time derivative of the local frame at 
fixed K using eqs. (20) and (26), the elements are 


m„, - X'Vijp - - «<*>**. 

(29a) 

m a 3 = - U^f, 

(296) 

(3) jt la) -(3) .(<*) _ « (3) tt (of) 

m 3a = e- Uije) - e i e, - 2e^ U, <J e j , 

(29c) 

(3 )rr (3) 

m 3 3 = e- Uijtj . 

(29 d) 

The rotation term R E is e^Uije^ if the polar axis is chosen as one of the eigen- 
vectors of the mean gradient matrix; its general expression is available in Cambon 
et al. (1985). The last equation relevant to our study is that which governs the 

pressure: 

(30) 

Without mean distortion, eqs. (28) and (30) correspond to a pure acoustic regime, 
where energy is exchanged between dilatational velocity and pressure at a frequency 
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ak. (The sonic speed a is easily reintroduced using the isentropic relation a 2 = 
yRT — jP/p-} On the other hand, the (exact) balance between the two last terms 
in Eq. (28) is the equivalent in physical space of the Poisson equation for p in the 
pure solenoidal flow. The solenoidal contribution to velocity is seen to be completely 
uncoupled from the dilatational field if m a3 is zero. This is valid for any irrotational 
compressing mean flow, but not for pure shear , as has been stressed by BMR. 
Finally, we note that the coupling of the solenoidal and dilatational fields is mediated 
by m 3a . This term is zero for spherical compression but must be considered for 
any anisotropic straining process (except for very specific wavevectors given by the 
particular deformation). An investigation of the timescales in (28) introduces the 
parameter R a (k) = ( t d ) for which Am is an averaged approximation in 

physical space. For very low values of this parameter, the incompressible limit is 
recovered, the dilatational mode tends to zero, and the sonic speed a approaches 
infinity; both kp/ p (which tends to the solenoidal solution to the Poisson equation) 
and its time derivative (which from (30) is observed to be proportional to o 2 p^) 
tend to finite non-zero values without inconsistency. At moderate R a (k ), a pseudo- 
acoustic regime is recovered, which deviates from the pure incompressible (itj >t = 0, 
P = Ps) case since the time variation of the m 3a term in (28) can be neglected 
and a WKB approximation can be used to predict the oscillating behavior of d 
(Sabel’nikov 1975 and DZ). (This oscillating behavior will be revisited in §5.) For 
large values of R a (k), the pressure term in Eq. (28) can be neglected compared 
to the other terms, and the “pressure released” regime is obtained. We note that 
use of the solenoidal Poisson equation to approximate the total pressure variance 
(i.e. setting p — p 9 ) and then using (30) to estimate the pressure-dilatation term, 
a method followed by DZ, can lead to some inconsistencies. If p = p 6 holds, the 
dilatational mode is directly given by a simplification of (28) (equating the first 
three terms to zero), and the solution is 


£ ( 3 ) ( M ) = j^ 3 \k,o). 


( 31 ) 


The potential inconsistency is that this solution for q 2 d is not necessarily the same 
as that found from Eq. (30). In (31), the dilatational part of the kinetic energy 
depends only upon its initial value; for a mean compression, q 2 d (t) = T{J{t))q 2 d ( 0), 
where f(J) depends on the type of compression. In contrast, the DZ method 
amounts to connecting both p and (p^ to the initial value of the solenoidal modes, 
^ o) (K,0), so that the dilatational part of the kinetic energy depends only on the 
solenoidal initial data: y^(t) = ^Dz(«J(0)(^ m ) 4 <?j(0) (where ^ ^*(J) and 

again depends upon the compression type). 

An approach which avoids this ambiguity and allows the classification of other 
relevant limits is available by introducing integrating factors into (28) and (30) so 
that 


V — J 


.^ 3 ) 


k 


and z = J 1 * 

pa 2 


(32) 
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FIGURE 1. Solenoidal and pressure-release regimes, 
satisfy the simpler equations: 

IMi . !> + e y = \z- ( 33 ) 

Vt 

+ g 2 z = a 2 z', (34) 

Vt 

where z’ = J^p./pa 2 = \{J / ka 2 )m 3a ft a \ The left-hand sides of both (33) and 

(34) are linked only to the dynamics of the solenoidal field and thus decoupled from 
the dilatational and pressure terms for irrotational mean deformations. 

We are now in a position to distinguish the different regimes implied by equations 

(33) and (34): 

I. The incompressible limit, with a 2 -+ oo, which corresponds to a vanishing value 
of all the time-derivatives in both equations; hence, z* and y-*0, and z = z* 
(i.e. p = pa) and y — 0 are consistent limits in this case. 

II. The acoustic regime, recovered when k 2 y > iz*. 

HI. The regime studied by Durbin & Zeman, where the pressure- dilatation correlation 
is given by the solenoidal pressure variance, which corresponds to k 2 y = iz* in 
(33) and z = z* in (34); these equalities hold only if the time- derivative of the 
solenoidal term (right-hand side of (33)) is much larger than the time-derivative 
of the dilatational term (first term on the left-hand side of (33)). 

IV. The pressure-released limit, corresponding to k 2 y < iz* in (33), which leads to 
the condition 2 

«Ji<( Am ) 4 fr 

(with L a lengthscale of the energy containing turbulence) required for the pres- 
sure-released regime to be valid. We mention in passing that if one assumes that 
the ratio A of the dilatational to solenoidal kinetic energy is proportional to M t , 
the above inequality suggests that an alternative to Am as the parameter that 
defines the pressure- released regime is the quantity AmM, 1/2 = rM t . In spite 
of this, the DNS results presented below indicate that the pressure-released limit 
seems to be adequately parameterized by Am alone. 
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Complete solutions of the system of line&r cqu&tions (28) end (30) are required for 
intermediate values of R a (k). The method described in Cambon (1982) of extending 
the Townsend (1976) approach can also be used for the general (homogeneous) case; 
an overview of the scheme follows, which we plan to follow in future work. A linear 
transfer function g, which generates the general solution of the system of equations 
as 

^(M) = 9ij (e { 3 \K/K< h R a (Ko),aoK 0 S) £< j) (K,0), 

is computed by solving (numerically in general) (20), (27), (28), and (30) for a 

set of arbitrary simple initial data (<^ ^ = 6 n , = 6 , 2 , etc ). (For convenience, 

the subscript can be taken to vary from 1 to 4 in order to represent the pressure 
fluctuation with = p/pa ) All the relevant one-point correlations can then be 
obtained by integrating over spectral space products of linear transfer functions and 
initial spectra. The initial spectra such as 

(^•>(P,0)^)(K,0)) = K - P) 

can be generated by invoking isotropy and assuming acoustic equilibrium (Sarkar 
et al 1989) and certain relationships between the dilatational spectra (Bataille, et 
al 1992). Note that the delta function S amounts to a factor <5 nm (A(0)/27r) 3 if 
discrete Fourier modes ( P = n,K = m) are chosen using periodic boxes of size A, 
as is done in DNS. For both the continuous and discrete case, the mean compression 
must be taken into account when computing correlations and integrating over wave 
space, using either 

tf(k-p) = *(K-P)J(f) 

or A (t)f A 3 (0) = J and <f 3 k = dk\dkidk^ = J _1 d 3 K. For the solenoidal field 
certain results can be obtained analytically, as is demonstrated below for the case 
of axial compression, since gfj depends only on the orientation of the wave vec- 
tor, and not on the modulus; integrations over wave space needed to derive the 
velocity correlations can thus be separated into the product of two one-dimensional 
integrals, one of which defines (independently of initial spectra shape) the initial 
kinetic energy. Evaluation of the non-solenoidal correlations is not as straightfor- 
ward since the components of the linear transfer matrix that involve the dilatation 
depend on both the direction e| 3) = k { /k of the wavevector k (as for the solenoidal 
case) and on its modulus. Accordingly, amplification coefficients like the functions 
T and Tbz mentioned above in general require numerical integration. This com- 
plication is a symptom of the wave number dependence of the sonic timescale in 
spectral space (a(O)A) 1 , symbolically shown in Figure la; since the deformation 
scale D 1 is the same for all wave numbers, above a critical value A"*, the sonic 
is the shorter of the two timescales. For a given energy spectrum with peak at Ao 
so that Am is characterized by A a (A"o), the rapid distortion behavior depends on 
K/Kq. The largest structures (A" < K 0 ) will, therefore, naturally tend toward the 
pressure- released extreme and the smallest (A" > K 0 ) toward the solenoidal limit. 


ullil.il i .11 



FIGURE 2. Contours of turbulent Mach number (a) before and (b) after axial 
compression (Case ClDW). 

When Kq falls well below 7\ *, the entire flow is within the pressure- released regime, 
and Am > 1; when K 0 > K*, the Am 0 limit is valid (see Figure lb). 

In the next section, the analysis is applied to the special case of axial compression, 
and DNS results are used to verify the relevance of Am as a critical parameter. 

3. RDT and DNS of axially compressed flow 

Both the RDT and DNS impose upon isotropic compressible turbulence the axial 
deformation that satisfies the homogeneity condition (15) so that the single nonzero 
mean velocity gradient component is 

U tl = — %- =DoJ~ l and F n = J. (35) 

For Do = D(0) < 0, this straining can be maintained for a finite time for as long as 
the flow volume is nonzero. Here we consider mean density ratios (equal to J~ ) 
that vary from 1 to 5; see Figure 2. Before describing the various DNS runs, which 
is done in §§3.2, in the next subsection we specify the RDT correlations relevant to 
one point modeling of the axial compression. 

S.l. Rapid distortion analysis for axial compression 
For the case of axial compression, the Craya-Herring-Cambon coordinates given 
in §2 reduce to e [ 3) = cos 0, e[ 2) = - sin (9, where 0 = (k, n) if the polar axis is chosen 
along the compression direction so that = 0 (see Cambon & Jacquin (1989) for 
other axisymmetric RDT applications). The RDT solutions for the solenoidal field 
are then 

£(D(k,0 = J£ (,) (K,0) and <^ 2 >(k,<) = j<? (2) (K, 0), 


(36) 
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with k\ = K i J~\ k 2 = K 2 , jfc 3 = tf 3 , cos* = k x /k = ZJ- l j{ 1 + C 2 £ 2 ) x / 2 , ifc 2 = 
1 + C 2 f 2 , C 2 = J” 1 - 1 and £ = K\/K. The double correlations are calculated 
using 

(u*u;) = and ux = — sin* + cos*. (37) 

Assuming isotropic initial data, both the solenoidal and pressure-released analytical 
RDT predictions can be obtained by integrating either over £ or directly in physical 
space, with the results being unaffected by the initial spectral shape. The axial 
compression correlations are tabulated below, using as super- or subscripts “s” and 
“p”, respectively, to denote the solenoidal and pressure-released limiting cases. 
Turbulent kinetic energy: 

= MJ), (38 a) 

= A P {J). (386) 

Compression-direction Reynold stress component: 


? 2 ( 0 ) 


KjfjKO 

9 . 2 ( 0 ) 




(39a) 


? 2 ( 0 ) 


B P (J). 


Compression-direction solenoidal-dilatational cross-correlation: 


(396) 


K“i>(<) = 0, 


(40a) 


( U M )(0 

9?(0) 


C P (J). 


Structure dimensionality tensor (Reynolds 1990; see also (47) and (48)): 


(406) 


DU(t) 

9»(0) 


= D.( J) = 


(41 a) 


Dnjt) 

? 2 ( 0 ) 


D P {J). 


Compression-direction component of the pressure-strain rate correlation: 


(416) 


i M) 
«JC o) 


DE,(J), 


(42 a) 


nn = e p {j ) = o. 


(426) 
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The pressure-dilatation term II is identically zero for both limits. In terms of the 
inverse compression ratio J with C 2 = J 2 — 1> we have for the J < 1 case. 


*-i 

( 2 tan 1 C\ 

( 1 + J C ) 

(43a) 

4 

2 + J 
ip ~ 3 

(436) 

B - J ' 2 

(1 _ C 2 — 1 tan -1 C \ 

(43c) 

' ~ 2C 2 

\2 + 2 C ) 


B - J — 
Bp ~ 3 

(43d) 

c P = ^( 

2 tan -1 C \ „ 

-1 + J 2 — - — - B, 

(43c) 

p 2 C 2 V 

C J 



tan -1 C\ (2 +J~ 2 \ 

(43/) 

E • - ( C1 

+ 3 + (C I -3)J- i ^!£). 

(43g) 


The solenoidal simplification functions turn out to be nearly linear in J 1 , whereas 
the pressure- released expressions are nearly parabolic. The quantities A„ and B„, 
previously derived by Ribner(1953), and the new expressions D, and E, are almost 
the same as those found for incompressible axisymmetric strain (of arbitrary his- 
tory), with, for example, the Reynolds stress tensor R*j = (u’uj) = J 2 ' 3 R*j(J), 
where R* is the RDT solution for the trace-free part of the mean deformation. 
Functions A p and B p are obtained by simply ignoring the pressure terms during 
the integration of the equations for the one-point correlations in the rapid axial 
compression limit, which are: 

\q 2 = -DRn + II, (44) 


Ru = — 2DRn +n u , 


(45) 


where Rij = ( UiUj ), and II = Tla/2 — (pn»,i)/P- 

For moderate compressibility, we find from the DNS results that the role of n u , 
which reduces the anisotropy (bn = Rn/q 2 — 1/3) in (44) and, therefore, indirectly 
reduces the production term in the kinetic energy equation, is more important 
than the direct role of II. For future reference, we now put forth some useful results 
concerning the anisotropy tensor bij = (uiUj)/q 2 — ^bij, and a recommended general 
decomposition for the Reynolds stress: 


Ru =«j(^+ »:{•’ + C’) + w«j> + <«f«s> +«s(^+ C) . 


(46) 
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where the solenoidal deviator b‘ } is split into a part reflecting the directional depen- 
dence (superscript e) and the “polarization” (subscript z) (Cambon et al. 1992). 
This decomposition can be recovered using the “structure tensors” introduced by 
Reynolds(1990) (see also BMR) via 

K = = j 2^|ie-(k,()rf J k = 26'<' ) ) (47a) 

c 'i = <a,hM = / + .*<•> - 6'J' 1 ) (476) 

Oii = few.;) = j 2^e'(k,i)<Ck = ,3 + i;]' 1 ) . (48) 

To show the equivalence of the two approaches we first note that the tensors in 
(47) and (48) are formally defined using the vector V>* or scalar ip potential func- 
tions according to the Helmholtz decomposition (18) in physical space. The three- 
dimensional spectra are then defined (the second equality in (47) and (48)), with e 9 , 
e d and being associated with respectively = 

I(^ 3 >*(^ 3 )), and (w*ujj). Finally, the directional-polarization anisotropy tensors 
are specified (the third equality) so that the two approaches are reconciled. 

In the notation of Reynolds (1990), D*j reflects the “dimensionality” of the 
solenoidal field, which is close to the directional anisotropy, whereas C’j is asso- 
ciated with the “componentality” of the turbulence. Only the dimensionality, or 
directional anisotropy, is needed for the dilatational velocity field since Df } = (ufu'j) 
(BMR). This reflects the single-component character of uf , in contrast to the two- 
component structure of the solenoidal field. We finally note that the above tensors 
are not all independent; for example, Cfj can be derived from the other tensors 
using (46)-(48), or equivalently the equation found by Reynolds (1992): 

(uW + Dr+Ctj^qfrj. 

It is hoped that the above general expressions will be useful in future attempts 
to model compressible flows. For the present, however, we narrow our approach as 
we use DNS results to test a few aspects of the rapid distortion analysis. 

S.2. Comparison to DNS of rapid axial compression 

The DNS results were obtained using a pseudo-spectral method to solve the 
compressible Navier-Stokes equations over a homogeneous domain in coordinates 
that move with the mean deformation (Rogallo 1981, BMR, Coleman k Mansour 
1991). As mentioned previously, the mean density ratio, J~ l = p(t)/p( 0) varies 
from 1 to 5 during the compression. The runs use for initial conditions compress- 
ible isotropic turbulence at various turbulent Mach numbers that have evolved from 
velocity fields, with finite dilatational components, that are in near acoustic equilib- 
rium; these initial fields are generated by running the code with no mean straining 
until they develop realistic triple-velocity correlations and dilatational energy for 
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FIGURE 3. Turbulent kinetic energy histories: , lower, Eq. (38a); , 

upper Eq. (386); , DNS with initial Am ranging from 0.3 (lower) to from 

3.0-8.0 (upper); , DNS Cases: lower, C1DJ (Af t , Am) t=0 = (0.03,0.3); upper, 

C1DV (0.1,7); middle, C1DW (0.3,1). 

the given M t . (Note that BMR have found that compressible isotropic turbulence 
strongly depends upon all the initial conditions for the dilatational field, not just 
Aft, which implies that had we begun the precomputation with, for example, a 
purely solenoidal field, the levels of dilatational energy in the developed flow might 
be significantly different than those found here.) The initial turbulent Mach num- 
ber for the runs varies from 0.03 to 0.44, the initial nondimensional compression 
speed r = \D\q 2 /e ranges from 50 to 800 (and |Z>|/(u>iWi) 1/2 from 2 to 88), and the 
initial values of Am = Af ( |D|/<u>,a>,) 1/2 fall between 0.26 and 7. A ratio of constant 
specific heats 7 = 5/3 and temperature dependent viscosity fx - T is ^sumed. 
All the runs used 96 3 grid points and were generated on the Intel Hypercube/i860 
at the NASA Ames Numerical Aerodynamic Simulation program. 

Results for the total (solenoidal and dilatational) turbulent kinetic energy will first 
be presented. In Figure 3, the DNS histories for (/>«.«, )/p are plotted against the 
mean density ratio J ~ 1 = p(t)/7>(0). (Because it is convenient in the code to solve 
for momentum rather than velocity, all of the DNS results presented approximate 
velocity correlations by using density weighted averages. We find for our purposes 
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1-5 2 Z* 3 3.5 4 4.5 5 


p(i)/p( 0) 

Figure 4. Solenoidal turbulent kinetic energy histories: , C1DJ (M t , 

Am) t=0 = (0.03,0.3); — ClDV (0.1,7); , ClDW (0.3,1); , Eq. (38a).’ 

that the uncertainty introduced by comparing the DNS Favre averages to the RDT 
Reynolds averages is unimportant.) These curves strongly support the validity of 
the RDT analysis presented above in that all the DNS results lie between the lower 
solenoidal (“4,”) and upper pressure-releases (“.4 p ”) RDT limits and that rate of 
energy amplification scales almost monotonically with the initial value of Am, which 
varies from 0.3 for the lower (dotted) curve to from 3 to 8 for the upper (dotted) 
curves. Three runs will be examined further, those represented by the dashed curves 
in Figure 3. Cases ClDJ, ClDV, and ClDW have initial Mt equal to 0.03, 0.1, and 
0.3, respectively, but the compression rates are such that the corresponding order 
for Am is 0.3, 7, and 1. At the end of the compression the ( M t ,Am ) values for 
ClDJ, ClDV, and ClDW are respectively (0.03,13), (0.2,79), and (0.4,6). 

Figure 4 confirms that Eq. (38a) is an excellent approximation for qj for the 
three cases considered and that the solenoidal field is, in fact, unaffected by the 
dilatational field, as predicted by the RDT. Both contributions to the kinetic energy 
are shown in Figure 5. We see that the dilatational energy is most important at the 
end of the compression, when the pressure-released regime dominates. The initial 
values of the dilatational- to-solenoidal energy ratio A 0 for the various runs is also 
apparent. 
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FIGURE 5. Solenoidal and dilatational turbulent kinetic energy histories. •, 

solenoidal; 0, dilatational; , C1DJ (M,,Am) M = (0.03,0.3); , C1DV 

(0.1,7); , C1DW (0.3,1). 

These results suggest the following model for the Mach number dependence of 
the kinetic energy behavior during a rapid axial compression. 

q](t) = A t (J)q]( 0) (49a) 

q](t) = A p (J)q 2 d (0) + (A+( J) - At(J)) q]{ 0 ), (496) 

where the “interpolation functions” A p and Af are assumed to vary monotonically 
with Am, increasing from zero to maxima of A p and A s , respectively. Similar agree- 
ment with DNS data is found for the other correlations given in (43). The results 
for (ujtt®) and {ufuf) are presented on Figure 6, where the DNS and RDT histories 
closely correspond. The slight overamplification of the DNS result compared to 
the analytical (u\u\)/q] = BJA, ratio becomes more pronounced with increasing 
Am. For the dilatational curves in Figure 6, (ufaj)/^’ we t ^ ie expected trend 
with Am, since they are closest to the analytical pressure-released expression (the 
“chain-dash” curve) when Am is largest. An analog to (49) is therefore proposed 
as a model for the dilatational Reynolds stress: 

(tifiif) _ o + B+(J) - BtjJ) 

*ld A p (J) Aq + A p ( J ) — At ( J) 


( 50 ) 
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FIGURE 6. Histories of the anisotropy of solenoidal and dilatational turbulent ki- 
netic energy: •, solenoidal; O , dilatational; , C1DJ (M t , Am) <=0 = (0.03,0.3); 

- — , ClDV (0.1,7); , ClDW (0.3,1); , Eqs. (38a), (39a); , Eqs. 

(38a, b), (39a, b), using A 0 = 0.22 from Run ClDV. 

where Aq is the initial ratio of the dilatational to solenoidal kinetic energy (which 
in practice might be neglected). The curves in Figure 6 suggest that the ratio 
(B+ — Bf )/(B p — B,) is smaller than the same ratio of “A” functions. 

Another anisotropy measure is investigated in Figure 7, where the structure ten- 
sors are presented. Recall that = {ufuf}. The fact that D, = B, in (41a) 
confirms that ^ b 3 tJ and = 1 b 3 } is a good approximation for ax- 

isymmetric strain, as suggested by studies of non-isotropic initial data under rapid 
rotation (Reynolds 1990, Cambon et al. 1992, Mansour e< al. 1991). Rapid rotation 
was shown to damp and, therefore, to reveal the initial anisotropy of as 
the asymptotic limit reached after several revolution times. In axisymmetric tur- 
bulence, Dn/q 2 can be interpreted as an angular coefficient cos 2 a, as implied by 
the integrands in (47) and (48), which reveals the conical structure of the spec- 
tral region that contains energy (around the symmetry axis). For example, a value 
of 1/3 for this coefficient suggests no angular dependence (directional isotropy), 
whereas a value between 0 and 1/3 suggests a relative concentration of spectral 
energy in the plane normal to the symmetry axis. Unfortunately, the situation is 
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FIGURE 7. Structure tensor histories: •, solenoidal; O , dilatational, , 

C1DJ (M,,Am )| s0 = (0.03,0.3); , ClDV (0.1,7); , ClDW (0.3,1); , 

Eqs. (38a), (41a); , Eqs. (38a, 6), (39a, 6), using A 0 = 0.22 from Run ClDV. 

more complex in the presence of a mean distortion, which causes a variation in di- 
rection of the time-dependent wavevector; in the pressure-released case, the angular 
distribution of spectral energy is unchanged with respect to (isotropic) initial data, 
but the wavevector tends to be aligned with the symmetry (compression) direction 
(see (36)) so that cos 2 a increases and tends to 1. On the other hand, in the pure 
solenoidal limit, the relative concentration of spectral energy in the plane normal to 
the compression direction opposes the tendency induced by the wavevector motion 
so that a slower (as compared to the pressure-released case), but still positive, net 
increase of the anisotropy is obtained. Note that the solenoidal ratio of j fq t 
given by the DNS is found to be slightly lower than the RDT analytical prediction. 

The cross- correlation (ujuf)/(ujtxj) is plotted in Figure 8 and compared to the 
RDT expression C p {J)/B a {J) from (39a) and (40b). The results suggest that for 
modeling purposes it might be advantageous to use an effective saturated vol- 
umetric ratio J + in place of J and define C+, an interpolating function for the 
cross-correlation, according to 

C P + (J) _ C„(J+) 

Bt(j) B.(j+y 


(51a) 
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p(t)/p(Q) 

FIGURE 8. Histories of the compression-direction component of the dilatational- 

solenoidal Reynolds stress correlation: , C1DJ (M t ,Am)/ =0 = (0.03,0.3); 

, C1DV (0.1,7); , ClDW (0.3,1); , Eqs. (39a), (406). 

and use the model 

= C p (J*)q 2 M (516) 

The parameter would tend toward the actual J in the pressure- released limit 
and approach unity in the solenoidal limit. The role of pressure will be discussed 
further in §4; for now, we observe in Figures 9 and 10 the dramatic increase of 
both pressure variance and pressure-dilatation terms caused by the compression. 
The amplification increases with the initial turbulent Mach number, which at first 
seems to conflict with the idea of a pressure released limit. The paradox disappears, 
however, if the pressure-dilatation term is no longer nondimensionalised by initial 
values (as is done in Figures 9 and 10), but rather scaled by a term proportional 
to the kinetic energy production. DNS results for Yl/Dq 2 are presented on Fig. 
11. The magnitude of this term is found to decrease with increasing Am for the 
three cases considered. This implies a non-monotonic variation with Am for this 
term (since it is identically zero in the solenoidal limit) with a maximum reached 
at low compressibility. It can be noticed that increasing values of Yl/Dq 2 are found 
at large J~ l for the intermediate Am case (ClDW), which we expect cannot be 
explained by RDT. This illustrates that the requirements for a compression to be 
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FIGURE 9. Pressure variance histories: , ClDJ (A/t,Am)t = o — (0.03,0.3), 

— ClDV (0.1,7); , ClDW (0.3,1). 

rapid enough for RDT to be valid are more difficult to meet when the flow is intrin- 
sically compressible, a fact also stressed by Zeman & Coleman (1992). The term 
Uu/Dq 2 linked to the compression- direction component of the pressure-stram rate 
correlations is shown in Figure 12. The solenoidal RDT expression, E t (J)/A,(J), 
from (38a) and (42a) is plotted and is found to give an upper limit to the DNS 
curves. These results suggest a monotonic decrease of Tlu/Dq 2 with increasing 
Am. Moreover, comparisons of the order of magnitude for both terms on Figures 
11 and 12 (noting the different scales of the two plots) show that the compression- 
direction component of the pressure- strain rate is dominant compared to its trace 
(pressure-dilatation term) in all cases. This confirms that the reduction of amplifi- 
cation of turbulent kinetic energy with respect to the pressure-released case (where 
only the “production” effects are present) is mainly due to Iln, through reduction 
of anisotropy, as in the pure solenoidal case. 


4. Towards a pressure-strain rate model 

Equations for Iln and II, valid for the rapid mean compression case, can be 
derived from eqs. (44) and (45), using eqs. (49) and (50) to model q 2 and (uju,). 
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FIGURE 10. Pressure-dilatation correlation histories: , C1DJ (Mt, Am) t =o = 

(0.03,0.3);- — , C1DV (0.1,7); , C1DW (0.3,1). 

The result is 

n„ = 

= (JHB. + B l+ - B+),’( 0) + J*B,qH 0)) 

= DE,,]( 0) + (. J\B + - ,J(0), (52) 

and 

n = (i(i. + A+ - A+) + D(B, + b; - Bt)) q]{ 0) 

= (U A t - A t) + d(b; - Bt)) q 2 A0)- (53) 

To obtain the above, the relations 

J~ 2 j t (J 2 B .) = E a - f t (J 2 B p ) = 0; \A , + DB. = 0; \A P + DB P = 0, (54) 

have also been used. 
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FIGURE 11. Rescaled pressure-dilatation correlation histories: , ClDJ 

(M u Am), =0 = (0.03,0.3); — C1DV (0.1,7); , C1DW (0.3,1). 

j.l Proposals for Second- Order Modeling 
Two simple ideas for constructing the Eq. (52) and (53) “interpolation” functions 
(denoted by a superscript “+”) are proposed: 

1. Using two functions of Am, passing monotonically from from 0 to 1 so that 
A+ - A+ = f\(A p - A t ) and B+ - B+ = f 2 {B p - B,)\ if the time- variation of 
the interpolation functions is neglected, this leads to the model 

n n = n;,(i-/ 2 ) (55a) 

n = (/* - fi)D(B p - B,)g 2 .( 0). (556) 

Note that f 2 > f\ is consistent with the sign of II found in the DNS results and 
with the interpretation of dilatational energy histories in Figure 5- 

2. Using a “saturated” volumetric ratio J + instead of the actual J in the evaluation 
of the interpolation functions with + superscripts, so that A + ( J) = A(J + ). The 
equation for J + would be 

j+ = U^J + -Cj^^(J + - 1 ), 


(55c) 
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FIGURE 12. Histories of correlation of pressure and compression-direction velocity 

gradient: , C1DJ (M t ,Am) (=0 = (0.03,0.3); , ClDV (0.1,7); , C1DW 

(0.3,1); , Eqs. (38a) and (42a). 


where Cj+ is a modeling constant. The sonic timescale-damping term would 
allow 7 + to saturate close to unity as the regime of the flow approaches the 
solenoidal limit. 

4-2 Testing a Second- Order Model 

From our analysis of the three DNS cases, we find that they are in the regime 
where the production and the rapid redistribution terms axe dominant. The contri- 
bution of the pressure-dilatation is about 10% of the production in the worst case. 
This leads us in our attempt to model the DNS results to adopt the first proposal 
of the previous subsection, and consider a linear (in 6,j) model for the solenoidal 
rapid part (see Shih et ai , 1990) of the redistribution term, taking 1 — / 2 (see 
Eq. (55a)) to be an exponential function of Am. The mean and Reynolds-stress 
equations reduce to: 




T, = 


Rij,t — RikU j y k RjkUi^k "I" $ij CXp( Am/C^m), 
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Figure 13. Reynolds stress history, compression-direction component: , 

ClDJ (M„Am) w = (0.03,0.3); — C1DV (0.1,7); , C1DW (0.3,1); ., 

Eq. (55d); no symbols, DNS. 

with 

= fa 4 {S*khj + Sj k bki - ISnnbmnSij) 

4- — (tokihj 4- Slkjbki), (55d) 

15 

where Am = y/SijSjiM t q 2 / e and we have set a = 2.523 (to be consistent with the 
model of Launder e< al. (1975)), and C&m = 40. The quantity f l»j = {Ui,j ~ Uj,i)/% 
is the mean rotation tensor. 

The development of the axial component of the Reynolds stress, R u , as predicted 
by the above model for the three cases considered is shown in figure 13. We find 
that this simple model, where the effects of the redistributive term diminish when 
Mt increases, compares well with the DNS data. The development of the turbulent 
kinetic energy (see Fig. 14) is also well reproduced, indicating that the effects of 
the pressure-dilatation are, in fact, weak compared to the production term. No 
attempt was made to optimize the constant C& m since the pressure-dilatation term 
was neglected. This term does play a role in the development of the flow, and C& m 
should be optimized in conjunction with a model for the pressure dilatation term. 
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FIGURE 14. Turbulent kinetic energy history: , ClDJ (M t ,Am)| =0 — (0.03, 

0.3); , C1DV (0.1,7); , ClDW (0.3,1); •, Eq. (55d); no symbols, DNS. 

5. Spherical compression and pure shear revisited 

5.1. Isotropic spherical compression 
In the presence of a mean spherical compression, with 

Uij = D6 ijt D=-^- = DoJ- 1 ' 3 F ij = J 1 ' 3 S tJ and k, = (56) 

the coupling term m 3cr in (27) and (28) has zero value. The evolution of the 
solenoidal kinetic energy is then easily found to be given by the amplification coef- 
ficient J“ 2 / 3 . For the dilatational field, eqs. (33) and (34) remain of interest now 
with their right-hand sides equal to zero (since p 3 ~ m 3a ). Even in the absence 
of the right-hand sides, a WKB analysis of the equations would not in general be 
appropriate because the timescale variation of a 2 and k 2 is not necessarily small 
with respect to the expected frequency ak of the oscillating system (depending on 
the value of Am). Blaisdell (1992, private communication) has recently found a 
solution free of WKB assumptions; its validity is restricted to values of 7 close to 
5/3, but a general analytical solution is possible (work in progress). If 7 = 5/3, 
k 2 and a 2 have the same J“ 2 / 3 time dependence, so simple solutions in terms of 
exp(±ia(0)k(t)t), where k(t) varies as in (56), can be obtained for y and 2. The 
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history of q 2 d can then be derived from the initial (uncompressed) dilatational field. 
Assuming that acoustic equilibrium holds for the initial conditions, one can write 

q 2 d (t) = J- 2,3 ql( 0 ), ( 5? ) 

which is the same variation found for the solenoidal energy. The acoustic equilibrium 
assumption is realistic but perhaps not necessary; the initial balance between kinetic 
dilatational energy and potential (pressure) energy allows oscillating terms to be 
dropped, but the same final result could also be reached after a certain elapsed time 
because of the damping behavior of integrals such as 


J F(K)exp(2\a(0)KJ 1 ^ 3 t)dK, 

where F is defined by the initial energy spectra and is nonzero only for flows out of 
acoustic equilibrium. (This behavior is similar to that found for the case of rapid 
rotation.) 

The above considerations show that an oscillating regime, more general than 
the pure acoustic one, is not inconsistent with the pressure-released limit and that 
the latter can be used to derive the same relationship (57) found via the acoustic 
equilibrium assumption. We thus find that the spherically compressed flow lends 
support to the general approach advocated in this paper. 

5.2. Pure plane shear 

The case of shear flow is particularly interesting because all the coupling terms, 
most notably m o3 and m 3o , are present. The crucial parameter in the absence of 
compression (J = 1) is the shear 5 = dUi/dx 2 . Under this deformation, eqs. (27), 
(33) and (34) become 

(p ] + Sj& 2) = S (( fc a + * 2 ) 1/2 ) y (58a) 


v_ 

vt 


i k ? {2) ) = ~ s i wTWfi ) * 2y 

•• 2,2 _C2 V f hith 1^1 W*) 

V + a 0 k y + S k2 y - S fc3 ) ' 


( 586 ) 

( 59 ) 


with 

Uij = S6n6 j2 , ki = Ki, k 2 = K 2 - K\St and k 3 = AY 

Here the polar axis is chosen to be in the gradient direction (n, = 6 l2 ). The two 
solenoidal an( i £( 2 ) components Eire very close to the set (o> 2 ,V 2 U2) used in 
linear stability analyses for decoupling, for example, the Orr-Sommerfeld equations 
for parallel flows (cf. Waleffe 1990). Even in the pure solenoidal case (where y = 
i^( 3 ) j k = 0), the present approach appears to be more tractable than are classic RDT 
approaches (Townsend 1976). Unlike for a purely irrotational mean deformation, 
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the presence of the new coupling terms (mediated by mi 3 = —Sk^/kk^, m 23 = 
Sk 3 /ki 3 in the above equations) makes the solenoidal field no longer independent 
of the dilatational component. In addition, this coupling introduces the new term 
S 2 (k 2 / k 2 )y in Eq. (59). The pressure- released approximation amounts to neglecting 
a 2 k 2 y compared both to this new term and to the solenoidal right-hand side in (59). 
The Am » 1 regime then implies that (in physical space), 

— ui(0) — Stu 2 ( 0), u 2 — u 2 (0) and u 3 = u 3 (0), 

and leads to quadratic amplification, with respect to St, of the kinetic energy (which 
is more rapid than the nearly linear amplification obtained by numerically integrat- 
ing the solenoidal RDT solution for + £< 2 )*^ 2 )) over k-space). Note 

that the inviscid solenoidal RDT solution for the vertical velocity component is 
given by T>(V 2 u 2 )/TH = 0 in physical space (corresponding to Eq. (58) with y = 0) 
so that a rapid decay of u 2 is found. On the other hand, u 2 is conserved in the 
pressure-released inviscid RDT limit. 

6. Recap and conclusions 

The objective of this analysis has been to develop a rapid distortion theory for 
homogeneous compressible turbulence at finite Mach number and then use that 
theory to explore some issues related to one-point compressible turbulence models. 
We have applied the analysis to the case of axial compression and found that DNS 
results confirm the RDT prediction of two distinct flow regimes, one for vanish- 
ingly small turbulent Mach number and the other for flows with negligible sonic 
and turbulent timescale variations compared to the mean distortion. The latter is 
referred to as the pressure-released regime (since the fluctuating pressure field can 
be neglected in the RDT for this limit) and is defined by large values of the product 
of Mt and the ratio of the turbulent to mean deformation timescales. For large 
values of this parameter, we find that the intrinsic compressibility of the turbulence 
is responsible for an increase in the growth rate of kinetic energy with increasing 
M t , an effect exactly opposite to that usually attributed to the compressibility. It 
would seem that the reduction in kinetic energy growth rate due to compressibility 
observed in previous compressible homogeneous DNS studies can be attributed to 
“slow” terms with nonlinear and dissipative origin, such as the “extra” dilatational 
dissipation associated by Zeman (1990) with eddy “shocklets.” In the future, we 
plan to perform systematic comparisons between compressible RDT (from numeri- 
cal solutions obtained by the method presented in §§2.4) and existing DNS to allow 
an accurate differentiation between the “rapid” and “slow” terms, which are found 
to have opposite trends with respect to the effect of compressibility on the kinetic 
energy growth rate. 

For the axial compression, analytic expressions for the correlations associated 
with one-point closures for both the solenoidal and pressure-released limits have 
been given. These expressions have been used to propose methods of interpolating 
between the two limiting RDT cases in models for the pressure-strain rate correla- 
tion, II jy and thus account for finite Mach number effects. 
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All the DNS results were obtained using the facilities of the Numerical Aerody- 
namic Simulation program. This work was partially supported by the Laboratoire 

de Mecanique des Fluides et d’Acoustique, Ecole Centrale de Lyon. 
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